#pragma warning(disable:4661)

#include "IcemVolume.h"
#include "Invert.h"
#include "boost/multi_array.hpp"
#include "boost/numeric/ublas/matrix.hpp"
#include "boost/numeric/ublas/vector.hpp"

#include "cuda.h"
#include "cuda_runtime.h"

#include "CudaIcemVolume.h"
#include <stdexcept>
#include <cstdint>
#include <algorithm>
#include <iostream>

using namespace boost;
//using namespace std;
using namespace boost::numeric::ublas;



template<class T>
IcemVolume<T>::~IcemVolume(){
	cudaFreeArray((cudaArray*)cudaVolume);
	cudaFree(colormap);
	for (int i=0;i<3;++i)
		cudaFree(cudaImage[i]);
}

template<class T>
void IcemVolume<T>::SetTalairachColormap(unsigned char  cm[3*2000]) {
	
	int ret = cudaMalloc(&colormap,3*2000*sizeof(unsigned char ));
	ret = cudaMemcpy(colormap,cm,3*2000*sizeof(unsigned char ),cudaMemcpyHostToDevice);
}

template<class T>
void IcemVolume<T>::SetColormap(unsigned char cm[3*256]) {
	int ret = cudaMemcpy(colormap,cm,3*256*sizeof(unsigned char ),cudaMemcpyHostToDevice);
}

template<class T>
void IcemVolume<T>::GetSliceSize(int dir, int* sizeX, int* sizeY){
	*sizeX = sliceSizeI[dir];
	*sizeY = sliceSizeJ[dir];
}


template<class T>
IcemVolume<T>::IcemVolume(int* size, T* data, double alpha, double brightness, double contrast, unsigned char colormap[3*256], float mat[4][4]):
	alpha(alpha),
	brightness(brightness),
	contrast(contrast),
	step(0.5)
{
	filterLimits.x=0;
    filterLimits.y=0;
	for (int i=0;i<3;++i) {
		shape[i] = size[i];
		sliceSizeI[i] = -1;
		sliceSizeJ[i] = -1;
	}
	int N = size[0]*size[1]*size[2];
	
	cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<T>();//8,0,0,0,cudaChannelFormatKindUnsigned);
	cudaExtent extent = make_cudaExtent(size[0]/**sizeof(T)*/,size[1],size[2]);
	int ret = cudaMalloc3DArray((cudaArray**)&cudaVolume,&channelDesc, extent);

	cudaMemcpy3DParms copyParams = {0};
	copyParams.srcPtr = make_cudaPitchedPtr((void*)data, size[0] * sizeof(T),size[0] ,size[1]); // was : sizeof(T)
	copyParams.dstArray = (cudaArray*)cudaVolume;
	copyParams.extent = extent;
	copyParams.kind = cudaMemcpyHostToDevice;
	
	ret = cudaMemcpy3D(&copyParams);
	//ret = cudaMemcpyToArray((cudaArray*)cudaVolume,0,0,(void*)data, sizeof(T)*N, cudaMemcpyHostToDevice);
	//ret = cudaMemcpy(cudaVolume,data,N*sizeof(T),cudaMemcpyHostToDevice);

	
	ret = cudaMalloc(&(this->colormap),3*256*sizeof(unsigned char));
	ret = cudaMemcpy(this->colormap,colormap,3*256*sizeof(unsigned char),cudaMemcpyHostToDevice);
    SetTransformationMatrix(mat);
  
}

template<class T>
void IcemVolume<T>::SetTransformationMatrix(float mat[4][4]){

   for (int i=0;i<4;++i) {
		for (int j=0;j<4;++j) {
			transformation[j][i] = mat[i][j];
		}
	}

	vector<double> pos2(4);
	matrix<double> invTrans(4,4);
	matrix<double> trans2(4,4);
	for (int i=0;i<4;++i){
		for(int j=0;j<4;++j){
			trans2(i,j)=transformation[i][j];
		}
	}
	InvertMatrix(trans2,invTrans);
	std::vector<vector<double> > corners(8);
	for (int k1=0;k1<2;++k1) {
		for (int k2=0;k2<2;++k2) {
			for (int k3=0;k3<2;++k3) {
				corners[k1*4+k2*2+k3] = vector<double>(4);
				pos2(0) = k1 * (double)shape[0];
				pos2(1) = k2 * (double)shape[1];
				pos2(2) = k3 * (double)shape[2];
				pos2(3) = 1;
				axpy_prod(trans2,pos2,corners[k1*4+k2*2+k3]);
			}
		}
	}
	vector<double> min(corners[0]);
	vector<double> max(corners[0]);
	for (auto it = corners.begin();it!=corners.end();++it){
		for (int j=0;j<it->size();++j) {
			min(j) = std::min(min(j),it->operator()(j));
			max(j) = std::max(max(j),it->operator()(j));
		}
	}
	vector<double> volSize(max - min);	

	for(size_t dir=0;dir<3;++dir) {	
		volumeSize[dir] = size_t(volSize(dir)/step);
		minCoordinates[dir] = min(dir);	
	}
}

template<class T>
void IcemVolume<T>::SetSliceSize(size_t size[3]){
	for (int dir=0;dir<3;++dir) {
		size_t newI,newJ;
		switch(dir) {
		case 0:
			newI = size[0];
			newJ = size[1];
			break;
		case 1:
			newI = size[0];
			newJ = size[2];
			break;
		case 2:
			newI = size[2];
			newJ = size[1];
			break;
		}
		if(newI!=sliceSizeI[dir] || newJ != sliceSizeJ[dir]) {
			cudaFree(cudaImage[dir]);
			int ret = cudaMalloc(&(cudaImage[dir]),newI * newJ * sizeof(short3));
			sliceSizeI[dir] = (int)newI;
			sliceSizeJ[dir] = (int)newJ;
		}
	}
}

template<class T>
void IcemVolume<T>::SlicePreparation(double pos[3], int dir, double minCoord_in[3]){
	
	matrix<double> invTrans(4,4);
	matrix<double> trans2(4,4);
	for (int i=0;i<4;++i){
		for(int j=0;j<4;++j){
			trans2(i,j)=transformation[i][j];
		}
	}
	InvertMatrix(trans2,invTrans);
	matrix<double> Tr(identity_matrix<double>(4));
	for (int i=0;i<3;++i)
		Tr(i,3) = minCoord_in[i];
	Tr(2-dir,3) = pos[2-dir];
	matrix<double> RX(identity_matrix<double>(4));
	RX(1,1) = 0;
	RX(2,2) = 0;
	RX(1,2) = -1;
	RX(2,1) = 1;
	matrix<float> RY(identity_matrix<float>(4));
	RY(0,0) = 0;
	RY(2,2) = 0;
	RY(0,2) = -1;
	RY(2,0) = 1;
	matrix<float> transform(4,4);
	matrix<float> temp(4,4);
	axpy_prod(invTrans,Tr,temp);
	if (dir == 0) {
		transform = temp;
	}
	else if (dir==1) {
		axpy_prod(temp,RX,transform);
	}
	else if (dir==2) {
		axpy_prod(temp,RY,transform);
	}
	for (int i=0;i<4;++i){
		for(int j=0;j<4;++j){
			sliceTransformation[i*4+j] = transform(i,j);
		}
	}
}

template<class T>
void IcemVolume<T>::SetFilterLimits(float lowerLimit, float upperLimit) {
 filterLimits.x = lowerLimit;
 filterLimits.y = upperLimit;
}


template class IcemVolume<uint8_t>;
template class IcemVolume<uint16_t>;